IMPLICIT MASS-MATRIX PENALIZATION OF HAMILTONIAN DYNAMICS WITH 
APPLICATION TO EXACT SAMPLING OF STIFF SYSTEMS 

PETR PLECHAC* AND MATHIAS ROUSSETt 

Abstract. An implicit mass-matrix penalization (IMMP) of Hamiltonian dynamics is proposed, and associated dynamical 
integrators, as well as sampling Monte-Carlo schemes, are analyzed for systems with multiple time scales. The penalization 
is based on an extended Hamiltonian with artificial constraints associated with some selected DOFs. The penalty parameters 
enable arbitrary tuning of timescales for the selected DOFs. The IMMP dynamics is shown to be an interpolation between the 
exact Hamiltonian dynamics and the dynamics with rigid constraints. This property translates in the associated numerical 
integrator into a tunable trade-off between stability and dynamical modification. Moreover, a penalty that vanishes with the 
time-step yields order two convergent schemes for the exact dynamics. Moreover, by construction, the resulting dynamics 
preserves the canonical equilibrium distribution in position variables, up to a computable geometric correcting potential, 
leading to Metropolis-like unbiased sampling algorithms. The algorithms can be implemented with a simple modification of 
standard geometric integrators with algebraic constraints imposed on the selected DOFs, and has no additional complexity in 
terms of enforcing the constraints and force evaluations. The properties of the IMMP method are demonstrated numerically 
on the Af-alkane model, showing that the time-step stability region of integrators and the sampling efficiency can be increased 
with a gain that grows with the size of the system. This feature is mathematically analyzed for a harmonic atomic chain 
model. When a large stiffness parameter is introduced, the IMMP method is shown to be asymptotically stable and to 
converge towards the heuristically expected Markovian effective dynamics on the slow manifold. 
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1. Introduction. This paper deals with numerical integration and sampling of Hamiltonian systems 
with multiple timescales. The main motivation is to develop numerical integration methods for the dy- 
namics which resolves certain selected fast degrees of freedom only " statistically" , and that can be also 
used to sample accurately the canonical equilibrium distribution. Furthermore, as an ultimate goal, one 
also seeks good approximation of dynamical behavior, at least at large temporal scales. 

Hamiltonian systems with multiple timescales typically appear in molecular dynamics (MD) simula- 
tions, which have become, with the aid of increasing computational power, a standard tool in many fields 
of physics, chemistry and biology. However, extending the simulations to physically relevant time-scales 
remains a major challenge for various large molecular systems. Due to the complexity of implicit methods, 
the time scales reachable by standard numerical methods are usually limited by the rapid oscillations of 
some particular degrees of freedom. Since the sampling dynamics has to be integrated for long times, the 
time-step restriction associated with fast oscillations/short time scales in molecular systems contributes to 
the high computational cost of such methods. However, the physical necessity of resolving the fast degrees 
of freedom in simulations is often ambiguous, and efficient treatment of the fast time scales has motivated 
new interest in developing numerical schemes for the integration of such stiff systems. 

The problem of integrating stiff forces is relevant both for the direct numerical simulation of the 
Hamiltonian dynamics, as well as for the less restrictive problem of designing a sampling scheme with 
respect to the canonical ensemble. Sampling from the canonical distribution can be achieved by Markov 
chain Monte Carlo (MCMC) algorithms based on a priori knowledge of possible moves combined with a 
Metropolis-Hastings acceptance/rejection corrector (a historical reference is [31]). For complex molecular 
systems, however, such global moves remain unknown in general, and sampling methods consists generically 
in using either a Hamiltonian dynamics integrator with a thermostat (e.g., a Langevin process), or its 
overdamped limit, a drifted random walk (Brownian dynamics) (see [7] for a review and references on 
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classical sampling methods). Brownian dynamics of systems with multiple time scales suffers from similar 
stability restrictions (see pUl ST] for some practical issues related to Brownian dynamics simulations in 
MD). 

Broadly speaking, one may start by recognizing two approaches to the numerical treatment of stiff 
systems: 

(i) Semi-implicit, multi-step integrators and their variants (e.g., the textbooks |361ll4j. or the review paper 

[lOj and references therein, [20j for Brownian dynamics), which attempt to resolve microscopic 
highly oscillatory dynamical behavior. 

(ii) Methods with direct constraints, where the highly oscillatory degrees of freedom are constrained to 

their equihbrium value (e.g., [36l [28l |40l [19] and references therein). 
In spite of their differences the common key feature of all these methods is to balance a trade-off between 
stability restrictions and implicit time-stepping form, or in other words, between the computational effort 
associated with small time steps, and the computational cost of solving implicit equations implied by the 
stiffness. 

Although constrained dynamics remove, in principle, the stiffness of the associated numerical scheme, 
it introduces new difficulties and numerical problems. As an approximation to the original dynamics it 
modifies important features of the system; most importantly, the original statistical distribution. The 
principal goal of the proposed method is to replace direct constraints by implicit mass-matrix penalization 
(IMMP), detailed in Section [J] which integrates fDOFs, but with a tunable mass penalty. The method 
designed in this way achieves the two goals: 

(i) from the dynamical point of view, the IMMP method amounts to an appropriate interpolation between 

exact dynamics and constrained dynamics considered in the second family of the methods men- 
tioned above. Moreover, a freely tunable trade-off between dynamical modification and stability 
is obtained. 

(ii) from the sampling point of view, the IMMP dynamics preserves the canonical equilibrium distribu- 

tion, up to a time step error and an easily computable geometric correcting potential. This leads 
to Metropolis Monte Carlo methods that sample exactly the canonical distribution. When us- 
ing Metropolis schemes, the forces arising from the geometric correcting potential need not be 
computed. 

The idea of adjusting mass tensors in order to slow down fast degrees of freedom goes back to ^3j. In this 
paper, the author proposes to modify the mass tensor with respect to the Hessian of the potential energy 
function in order to confine the frequency spectrum to low frequencies only. Two natural drawbacks of this 
procedure arise from the costly computation of the second-order derivatives of the potential, and from the 
bias introduced when the adjusted mass-tensor is adapted during the dynamics. Such an approach seems 
inevitable when the fDOFs are unknown, but in many cases, the fast degrees of freedom are explicitly given 
by the structure of the system (e.g., co-valent and angle bonds in molecular chains). To our knowledge, 
mass tensor modification have been used in practical MD simulations by increasing the mass of some well- 
chosen (e.g., light) atoms [26l[29]. The aim of this paper is to propose a more systematic mass-tensor 
modification strategy. 

The proposed method relies on the assumption that the system Hamiltonian is separable with quadratic 
kinetic energy 

H(p,q)^^p'^M-'p + Viq), (1.1) 

and that the "fast" degrees of freedom {£,i,.-,£,n) are explicitly defined, smooth functions of the system 
position 

q=iqu...,qd)^ {^l{q),-,Uq)) ■ (1-2) 

We emphasize that the knowledge of "fast forces" is not required, and the variables ^ can be chosen 
arbitrarily. If the fDOFs are not identified the method retains its approximation properties while not 
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performing efficiently. The fDOFs are penalized with a mass-tensor modification given by 



M,(q) - Af + z.2v,eM.V^C, 



(1.3) 



where v denotes the penalty intensity, and a "virtual" mass matrix associated with the fDOFs. The 
modification does not impact motions orthogonal to the fDOFs. The position dependence of the mass- 
penalization introduces a geometric bias. This bias is corrected by introducing an effective potential 



which will turn out to be a z^^^-perturbation of the usual Fixman corrector (see |19| ) associated with the 
sub-manifold defined by constraining the fDOFs ^. The key point is then to use an implicit representation 
of the mass penalty with the aid of the extended Hamiltonian 



The auxiliary degrees of freedom z are endowed with the "virtual" mass- matrix M^. The constraints (Ci/) 
are applied in order to identify the auxiliary variables and the fDOFs with a coupling intensity tuned 
by V. The typical time scale of the fDOFs is thus enforced by the penalty v. The system is coupled to 
a thermostat through a Langevin equation (|2.3p . which yields a stochastically perturbed dynamics that 
samples the equilibrium canonical distribution. We then obtain the following desirable properties: 

1. The associated canonical equilibrium distribution in position is independent of the penalty v. 

2. The limit of vanishing penalization {y — 0) is the original full dynamics, enabling the construction 
of dynamically consistent numerical schemes. 

3. The limit of infinite penalization is a standard effective constrained dynamics on the "slow" man- 
ifold associated with stiff constraints on ^. 

4. Numerical integrators can be obtained through a simple modification of standard integrators for 
effective dynamics with constraints yielding equivalent computational complexity. 

The dynamics associated with the IMMP Hamiltonian (|1.5p is detailed in p.Sp , see Section [31 The 
numerical discretization (using a leapfrog/Verlet splitting with constraints, usually called "RATTLE" for 
fully constrained dynamics) is given by (|4.ip . When considering sampling, the time-step error of the 
numerical flow can be corrected with a Metropolis step (the so-called Generalized Hybrid Monte-Carlo 
method, see references in Section 14]) to obtain exact sampling. When this correction is introduced, the 
gradient of the Fixman potential (II. 4p need not be computed. These numerical aspects of the method are 
presented in Section 3] By including a penalty, the proposed method modifies the original Hamiltonian. 
However, the mass penalty can be also thought of as depending on the time step v = v{5t) leading to order 
two consistent schemes. 

In Section [6l we introduce a small stiffness parameter e encoding the fastest DOFs, and show that 
the penalty intensity can be scaled as v — v/e in order to obtain asymptotically stable dynamics in the 
limit e — + 0. We prove that the dynamics converge towards the expected Markovian effective dynamics on 
the slow manifold. We also present analysis of the corresponding asymptotic preserving properties of the 
proposed scheme. 

High-dimensional systems usually contain a large variety of timescales, and are therefore challenging 
test cases. The A^-alkane model is numerically studied in Section [5] and systematically compared to Verlet 
scheme and constrained integration, with separate studies for dynamical and sampling issues. In the 
case of butane, bond angles are penalized. Dynamical interpolation between exact dynamics and rigidly 
constrained dynamics is demonstrated, with the associated gain in the time step stability. On the other 
hand, exact sampling with possible gain in the mixing time for Metropolized sampling methods is analyzed. 



Va^M) = :T75l^(det {MM)) 



(1.4) 




(a) 



(1.5) 
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Furthermore, for the A^-alkane model with large iV, torsion angles are penalized, a large mass-penalty of 
order 0{N), i.e., v = DN, where N is the system size is considered. For dynamics, it induces a gain in 
the time step stability region that grows with TV, while numerical evidence is given that some macroscopic 
timescales, e.g., low frequencies of the chain length dynamics, remain of order 0(1). For sampling, it 
induces a similar gain for mixing time in terms of iteration steps, and measured with autocorrelation of 
the chain length evolution. Rigorous proofs with explicit scalings of this behavior are provided for the case 
of a linear atomic chain with quadratic (harmonic) interactions in Section [7l and consistence of the IMMP 
macroscopic dynamics towards a stochastic wave equation when the re-scaled penalty vanishes D Q is 
demonstrated. 

Acknowledgments: The research of M.R. was partially supported by the EPSRC grant GR/S70883/01 
while he was visiting Mathematics Institute, University of Warwick. The research of P.P. was partially 
supported by the National Science Foundation under the grant NSF-DMS-0813893 and by the Ofhce of 
Advanced Scientific Computing Research, U.S. Department of Energy; the work was partly done at the 
ORNL, which is managed by UT-BattcUe, LLC under Contract No. DE-AC05-00OR22725. 

2. Langevin processes and sampling of canonical distribution. We consider a Hamiltonian 
system in the phase-space x with the Hamiltonian H in the form 

H{p,q)^]^pTM-^p + V{q), (2.1) 

We use generic matrix notation, for instance, the Euclidean scalar product of two vectors pi,P2 G 

is denoted by pjp2, and the gradients of mappings from R'' to M" with respect to standard bases are 

represented by matrices 



(V[e)», - {VgOJ^ - 1^ , i - 1, . . . , n , i = 1, 



When the system is thermostatted, i.e., kept at the constant temperature, the long time distribution 
of the system in the phase-space is given by the canonical equilibrium measure at the inverse temperature 
/3 (also called the NVT distribution) given by 

Hidpdq)^ ^e~'^"<-P''^Updq, Z=[ e-'^"^P''^Up dq , (2.2) 

with the normalization constant Z < oo. The standard dynamics used to model thermostatted systems 
are given by Langevin processes. 

Definition 2.1 (Langevin process). A Langevin process at the inverse temperature /? with the Hamil- 
tonian H{q,p), {p,q) (z R'^ X Mf^, the d x d dissipation matrix j, and the fluctuation matrix a is given by 
the stochastic differential equations 



WgH --fq + aW, 



(2.3) 



where W is a standard white noise (Wiener process), and a € x M*^ satisfies the fluctuation-dissipation 
identity 

For any 7, the process is reversible with respect to the stationary canonical distribution 112. S^) . Furthermore, 
if 7 is strictly positive definite, the process is ergodic. 
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Throughout the paper, stochastic integrands have finite variation thus the stochastic integration (e.g., 
Ito or Stratonovitch) need not be specified. Furthermore, the usual global Lipschitz conditions (see [35]) 
on H and ^ are assumed, ensuring well-posedness of the considered stochastic differential equations. The 
analysis presented in the paper can be generalized to a position dependent dissipation matrix 7 = 7(9). 

The mapping ^ : R"^ ^ R" , defines n < d degrees of freedom, given by smooth functions taking values 
in a neighborhood of 0. We assume that the mapping ^ is regular (i.e., with a non-degenerate Jacobian) 
in an open (5-neighborhood Os = {q \ \\£,{q)\\ < S} of ^""'^(0), hence defining a smooth sub-manifold of R'^ 
denoted Mz = for z in a neighborhood of the origin. The dependence of the potential V with 

respect to the degrees of freedom ^ is expected to be "stiff" in the second variable. In Section [S] we will 
introduce the stiffness parameter e. In that section we shall assume that such parameter dependence can 
be explicitly identified, and that the potential energy V can be written in the form 



where the function U : M.'^ x M" R satisfies the coercivity condition lim^^+oo U{q,z) — +00. The fast 
degrees of freedom ^ of states at a given energy then remain in a closed neighborhood of the origin as the 
stiffness parameter e ^ 0. In this limit the system is confined to the sub-manifold Mq which is usually 
called the "slow manifold" . 

3. The implicit mass-matrix penalization method. In this section we focus on properties of 
the IMMP method. The multiscale structure of the potential V need not be known in order to apply 
the method. Thus, in this section, we consider the potential V in the form where we do not impose the 
structural assumption (|2.4p on the potential function y : R"^ ^ R. 

3.1. Description of the method. The new, penalized mass-matrix of the system is the position 
dependent tensor defined in (jl.3|) . The associated modified impulses are denoted 



When v becomes large, the velocities are bound to remain tangent to the manifolds {q\£,{q) = z}, and 
orthogonal motions are arbitrarily slown down. Conversely, when v — 0, one recovers the original highly 
oscillatory system. Since the modification in M^, depends on the position q, new geometry is introduced 
and an additional correction (ll.4|l in the potential energy is required in order to preserve original statistics 
in the position variable. This correction is in fact close to the standard Fixman corrector for v large (see 
(|6.5p ). Defining G{q) as the n x n Gram matrix associated with the fast degrees of freedom 




(2.4) 



(3.1) 



G{q)=V^^M-'V,i, 



(3.2) 



one has the following property of the correcting potential. 
Proposition 3.1. Up to an additive constant, we have 




(3.3) 



and thus (up to additive constants) 



lim Vfix.iy = Vfix — TT^li^dct {G{q)) , and 



lini Vfix.z. = . 



Proof. Using the identity for a non-diagonal matrix J of dimension ni x 712: 



det (Id „i + J J^) ^ det (Id + J'^ J) , 
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one observes 

det (M^) = dot (M) det (v^M^) det (G + ^M"^) 

from which the expression for the corrected Fixman potential follows. □ 
The associated modified Hamiltonian is then given by 

HAp., q) = IpIM,-'p, + V{q) + Va^Al) , (3-4) 

and Hq = H is the original Hamiltonian (jl.ip . 

Sampling such a system can be done using the standard Langevin stochastic perturbation as detailed 
in Definition 12.11 However, the direct discretization of the equation of motion given by H^, (e.g., by an 
explicit scheme) is bound to be unstable from non-linear instabilities when the fast degrees of freedom are 
not affine functions. In order to construct stable schemes one may rather use an implicit formulation of the 
Hamiltonian p.4p . in conjunction with a solver which enforces the constraints. To obtain such a formulation 
we extend the state space with n new variables and associated moments (p^^ , The 

auxiliary mass-matrix for the new degrees of freedom is then given by ■ The new extended Hamiltonian 
of the system i?iMMP, defined by (jl.Sp . is now defined in x where n position constraints denoted 

by (Cy) are included. This construction implies n hidden constraints on momenta. The equivalence of the 
two Hamiltonians (|3.4p and (|1.5p formulations is stated as a simple separate lemma. 

Lemma 3.2. The equations of motion associated with the penalized mass-matrix Hamiltonian (|3.4[) or 
the extended Hamiltonian with constraints (|1.5p are identical. 

Proof. The Lagrangian associated with 7?immp is given by 

LiMMp{q,z,q,z) = i(7^Afg+ ^z'^M^J - V{q) - Ffix,z.(g) , 

and includes hidden constraints on velocities i = i^^^f^q implied by the constraints (0,^) on position 
variables. Replacing z and z in £immp by their expressions as functions of q and q, one obtains the 
Lagrangian associated with H^. □ 

The stochastically perturbed equations of motion of the Langevin type associated with ()1.5|) define the 
dynamics with implicit mass-matrix penalization. 

Definition 3.3 (IMMP). The implicit Langevin process associated with Hamiltonian iJiMMP md 
constraints (C^) is defined by the following equations of motion 

■ q = M~ V 
z = M~^p, 

p = -%T/(g) - V,Ffi,,,((z) --iq + aW- A 
A 

Pz = -IzZ + a^Wz + - 

V 

The process W (resp. Wz ) is a standard multi- dimensional white noise, 7 (resp. ^z) a dx d (resp. nx n) 
non-negative symmetric dissipation matrix, a (resp. cTz) is the fluctuation matrix satisfying aa"^ — ^7 

(resp. <Jz<^'^ — ^Iz)- The processes A £ M" are Lagrange multipliers associated with the constraints {C^ 
and adapted with the white noise. 

This process is naturally equivalent to the explicit mass-penalized Langevin process in M'^ x M'' as- 
sociated with M,y. Moreover, when the penalization vanishes {ly 0), the evolution law of the process 
{Pt,qt}t>Q or {{p^)t, qt}t>o converges towards the original dynamics. 
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Proposition 3.4. The stochastic process with constraints (|3.5p is well-posed and equivalent to the 
Langevin diffusion (see Definition \2.1]) . with the mass-penalized Hamiltonian (j3.4p . and the 

dissipation matrix given by 

Furthermore, the process is reversible and ergodic with respect to the canonical distribution (|3.8p (with 
marginal in position variables given by the original potential, i.e., up to the normalization, e~^^'^'^^dq) . 
Proof. Imposing the constraints implies V^^il/^-'^p = ^AI~^pz. Thus by the definition of p^ we have 

Pi. = p + vVq£, Pz . 

Since the position process {qt}t>o is of finite variation, a short computation shows that for each coordinate 
i = 1, ... d 



Furthermore, 



and thus 



Pl ^ f + ud.^ip, + v'q' ^13,^0 Pz ■ (3.6) 



P.^p + v^qip, - V, Qp^M- V.^ • 

Substituting the expressions for p and Pz from (|3.5p into p.6p we obtain 

= -^P^VqM- V - ^qV{q) - V<,Ffi^,^(g) - 79 - i/V,C7.2 + ctW^ + J^%ecr^W-^ , (3.7) 
which yields the result. □ 

3.2. Exact sampling in position variables. By construction, statistics of positions q of the mass 
penalized Hamiltonian are independent of the penalization, leading to the exact canonical statistics in 
position variables. 

Proposition 3.5 (Exact statistics). The canonical distribution associated with the mass-penalized 
Hamiltonian p.4p is given by 

H^idp, dq) = ^e-P^-^P-'^Up, dq . (3.8) 
Its marginal probability distribution in q is 

1 ( ,~^H.,,..)dp. ^''""''''^ 



which is the original canonical distribution (|2.2p in the position variables, and is independent of the mass 
penalization parameter v. 

Proof. The normalization of Gaussian integrals in the py variables yields 

d/2 



which is cancelled out by the Fixman corrector Vfix,i/ and the result follows. □ 
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3.3. Interpolation between exact and constrained dynamics. In this section, the IMMP dy- 
namics is shown to be an interpolation between exact dynamics =^ 0), and constrained dynamics 

(jy = +oo). 

Proposition 3.6 (Small penalty). When v ^ the evolution law of the processes {pt, Qt}t>o or 
{(p^)t, (7(}(>o defined by the implicit equations (|3.5p converges (in the sense of probability distributions on 
continuous paths endowed with the uniform convergence) towards the process solving the original Langevin 
dynamics (|2.3p . 

Proof. The stochastic differential equation defined by q = M^^p^, and p.7p has smooth coefficients 
which depend on in a continuous fashion [v i— > AI^, and i' i— > Vsx.i^ are continuous). Standard results on 
weak convergence ([16]) of stochastic processes imply the result as stated. □ 

When the mass penalty tends to infinity, the IMMP process converges to a constrained process on the 
manifold Mzt^o = U I ^(q) = ^t=o}- 

Proposition 3.7 (Large penalty). Consider a family of initial conditions indexed by v and satisfying 

sup {Ult=o) - zt=o)\ < +00 , 



and assume that the Gram matrix G is invertible in a neighborhood of M.zt=o- Then when v +oo the 
IMMP Langevin stochastic process (|3.5p converges in distribution towards the decoupled limiting processes 
with constraints 



q^M-^p, 

P = -VqV - VgVfix 

C(g) = zt^o , 
z 

Pz = ~lzZ + (JzWz ■ 



jq + aW - Vg^A , 



(C) 



(3.9) 



where {Xt}t>o are adapted stochastic processes defining the Lagrange multipliers associated with the con- 
straints (C). 

Furthermore, the process {qt,Pt\t>o defines an effective dynamics with constraints (see also Defini- 
tion \6.S\) on the sub-manifold Mzt=o- reversible with respect to its stationary canonical distribution 

given, up to the normalization, by the "stiff" Boltzmann distribution 



-[3(H(p,q)+VtUq)) 



^T'M^^^Jdpdq) 



with the q-marginal e"''^'*) 5^(g)=o('^'Z)- When 7 and 7^ are strictly positive definite the process is ergodic. 

Proof. By a simple translation, it is sufficient to show the proposition for zt^Q — 0. Satisfying 
the constraint (C^) in (|3.5p implies a hidden constraint in the momentum space, Vg^M^^p = ^M~^pz. 
Differentiating this expression with respect to time and replacing the result in (|3.5p yields an explicit 
formula for the Lagrange multipliers 



A=(G+^M-i)- 



Hess (0 (A/- V, M-'p) + V,^M-\f, - ^M'^fz 



(3.10) 



with forces (/g, fz) 



-\IqV - VgVfix,^ - 7il/- V + aW , 
-IzM-^Pz + crzWz, 



and the Hessian Hess (^) of the mapping ^ acting on the velocities M ^p. This calculation shows that 
is in fact a standard stochastic differential equation with smooth coefficients, and thus has a unique strong 
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solution. The coefficients of tfiese stocliastic differential equations are continuous in the limit ^ — > 0, at 
least in a (5-neighborhood of A^o in which G is invertible. The formally computed limiting process is given 
by (|3.9p with the Lagrange multipliers solving 

A = (Hess (0 (M- V, M" V) + V,) ■ 

By construction, this limiting process satisfies the constraint ^(q) = 0. Its coefficients are Lipschitz and 
the process is well-posed. As a result of those properties, the rigorous proof of weak convergence follows 
classical arguments, see jl6j . that are divided into three steps 

(i) We truncate the process (|3.5p to a compact neighborhood oi Mq. 

(ii) The continuity of the Markov generator with respect to i implies tightness for the associated ^- 

sequences of truncated processes with the limit being uniquely defined by (|3.9p . 

(iii) The limiting process remains on A^o, which implies weak convergence of the sequence without trun- 

cation. 

The process p.9p is thus a Langevin process with constraints, exhibiting reversibility properties with respect 
to the associated Boltzmann canonical measure and is ergodic when 7 is strictly positive definite (see the 
summary in Appendix [B| . Note that the g-marginal is geometrically corrected by the Fixman potential 
term. □ 

We conclude this section by discussing some consequences for numerical computations. 

Remark 3.8. Proposition l3 .61 and I3.7l implv that the IMMP scheme is a tunable interpolation between 
the exact stochastic dynamics (|2.3|) . and the stochastic dynamics with constraints p.9p . If one prefers to 
interpolate with rigidly constrained dynamics (i.e., without the Fixman correction, see also Section [5] for 
a detailed discussion on "stiff", as opposed to "rigid", constrained dynamics), one removes the Fixman 
correction in the force evaluation. 

4. Numerical integration. The key ingredient for achieving efficient numerical simulation is to use 
an integrator that enforces the constraints associated with the implicit formulation of the mass penalized 
dynamics p.Sp . The implicit structure of (|3.5p leads to numerical schemes that are potentially asymp- 
totically stable in stiff cases (Section [B]). On the other hand, when the penalization v vanishes with the 
time-step the scheme becomes consistent with respect to the original exact dynamics (|2.3p . One may then 
consider the mass-penalization introduced here as a special method of pre-conditioning for a stiff ODE 
system with an "implicit" , in the time evolution sense, structure. Here, the "implicit" structure amounts 
to solving the imposed constraints ^(g) = z/v 'm (|3.5p . 

It lies outside the scope of this paper to review standard numerical methods for constrained mechanical 
systems, we refer to [M] as a classical textbook, and to the series ( [40] [25l [TSl [2T1 [8]) as a sample of 
works on practical developments of numerical methods. The IMMP method is presented with the classical 
leapfrog/Verlet scheme that enforces constraints, usually called RATTLE in (|4.ip . It can be implemented 
by a simple modification of standard schemes constraining fDOFs. The scheme is second order, reversible 
and symplectic. This choice is largely a presentation matter, for practical purposes one can refer to one's 
favorite numerical integrator for Hamiltonian systems with or without stochastic perturbations. 

For accurate sampling of the equilibrium distribution, one can also add a Metropolis acceptance, re- 
jection time-step corrector at each time step of the deterministic integrator. If the underlying integrator 
is reversible and preserves the phase-space measure, this extension leads to a scheme which exactly pre- 
serves canonical distributions. The Metropolis correction is used in Hybrid Monte-Carlo (HMC) methods, 
which are sampling algorithms relying on the underlying dynamics of the system to generate moves in 
the configuration space, which are accepted or rejected according to the Metropolis rule. The Metropolis 
acceptance/rejection step can be used at each integration time step of the Langevin process, usually re- 
ferred to as Generalized Hybrid Monte Carlo (GHMC) introduced in [22]. However, the necessity of the 
momentum flip when a rejection occurs destroys the dynamical features of the Langevin process when the 
rate of rejection does not vanish, and makes the latter rather behave similarly to an overdamped dynam- 
ics. Note also that for usual schemes integrating the Hamiltonian dynamics, the average acceptance ratio 
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often decreases when the dimension of the system increases ([24j). Many improvements and modifications 
( [TTJ [551 [21 [23] ) of the HMC algorithm have been developed since its introduction in [T^] for simulations 
applied to quantum statistical field theories. Subsequently it has been also employed to a wide range of 
simulations in macromolecular systems, e.g., [9l 138). 

We implement numerical discretization of the Langevin process with constraints p.5|) obtained by 
splitting the Hamiltonian part and the Gaussian fluctuation/dissipation perturbation. Note that by using 
a Metropolis acceptance/rejection rule (HMC), or simply by weighting statistical averages, the forces 
associated with the Fixman corrector need not be computed when one is interested in sampling only. 



4.1. The IMMP integrator. We recall that we consider the IMMP dynamics (|3.5p . which consists 
of the following elements 

1. the Hamiltonian i/iMMP defined in ()1.5p . which defines the deterministic dynamics of the IMMP. 

2. the dissipation matrix diag(7,72), and the inverse temperature /?, which defines the features of 
the stochastic thermostat. 

Scheme 4.1 (Dynamical integrator). 

Step 1: Integrate the Hamiltonian part with: 





St, 


Pn+1/2 

< 


-Pn~ji 


^ Pn+1/2 


= Pn + ->^ 
V 


j qn+1 = 


qn + 6tM^ 


[ Z,i+1 — 


Zn + 5tM~ 




•j _ Zn+l 
V 




6 


Pn+1 = 


Pn+1/2 



Pn+1 ^Pn+1 + 



Vj^(q„+l)M-Vn+l = -M-Yn+1 



(C1/2) 



(4.1) 



iCi). 



Step 2: Integrate if necessary the Gaussian fluctuation/dissipation part with a mid-point Euler scheme with 
constraints (see Appendix\D^. 

Here again, one can remove the Fixman correction forces VqVfix.y in force evaluation in order to 
interpolate with usual rigid constraints dynamics. 

Note a useful variant of the integrator, which occurs when the potential dependence with respect to the 
penalized variables is known explicitly V{q) = U{q,£,{q)). In such a case, the expressions (|4.ip is replaced 
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by 



Pn - 77-^2U{qn) + -K+l/2 



Pn+1/2 




V 



iCi/2) (4.2) 



V,Vfix,,.)((7,i+i) - Vq^((7„+i)A„+i 
^V2J7(g„+i) + -A„+i 
W;f^{qn+i)M-^Pn+i = -Af-ip^+i {Ci 



where in the above, Vi and V2 denote the derivatives with respect two the first and the second variables, 
respectively. This variant consists in applying a part of the force to the auxiliary variables instead of 
directly to the system. For penalized degrees of freedom that move far from their equilibrium value, e.g., 
torsion angles, see Section [5l this may increase the stability of the algorithm enforcing the constraints. In 
the same way for slow/fast systems (see Section |6^ . the scheme (|4.2p will lead to asymptotic stability in 
the large stiffness limit. 

4.2. Exact sampling and Monte-Carlo scheme. To obtain exact sampling by correcting time 
step errors and, if necessary, the geometric bias, a Metropolis acceptance/rejection is added. A domain 
of phase-space Dst C IR^"^ x R^" where the integrator (|4.ip with constraints has a unique solution is also 
considered. In practice this set is simply the set of configurations for which the algorithm used to enforce 
constraints in (j4.ip (e.g., a Newton algorithm) converges in a given number of steps. In practice, the 
(equilibrium) probability for the system of lying outside Dst goes to zero exponentially fast with the time 
step, and the latter is taken such that this probability is negligible. 

Scheme 4.2 (Leapfrog/Verlet algorithm with Metropolis correction for Langevin IMMP (|3.5p ). 
Step 1: Compute (qn+i, Zn+i,Pn+i,Pn+i) '^^^^ tf'-^ integrator (14. ip . and set 

A-ff„+i = HiMMp{qn+l, Zn+l,Pn+l,Pn+l) ^ HiMMp{qn, Zn,Pn,Pn) ■ 

If {qn+1, Zn+i,Pn+i,Pn+i) does not belong to Dst, set AH„+i = +00. 
Step 2: Accept the step with the probability min(l, 6"'^'^^"+'^), otherwise reject, flip momenta, and set 

{qn+l,Zn+l,Pn+l,Pn+l) (?« , 2„ , -p„, -p^) . 

Step 3: Integrate the Gaussian fluctuation/dissipation part with a mid-point Euler scheme (for details see 
Appendix\D^ . 

Remark 4.1. A useful variant, when using HMC strategies, consists in modifying the Hamiltonian 
-ffiMMP (jl.Sp in the integrator (importance sampling) by neglecting the Fixman corrector 

HiMMp{p,Pz,q, z) = ]-p^Ar^p + ipf + y{q) 

2 2 (43) 

?(?) = l (a), 

where is a potential that may be chosen arbitrarily (typically V — V). Indeed, only the underlying phase 
space structure, which does not depend on V , is necessary in HMC methods. The correct potential V+Vfix,„ 
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has to be used in the Metropohs step only, or alternatively by weighting ensemble averages, in order to 
ensure exact canonical sampling. Thus potentially costly evaluations of the gradient of the Fixman corrector 
Vfix.i' can be avoided. This numerically constructed Markov chain preserves the canonical distribution. 

Proposition 4.2 (Exact sampling). The Monte-Carlo algorithm generated from (j3.5p as described in 
Scheme \4-^ generates a Markov chain that leaves invariant the canonical distribution (j3.8p ( conditioned in 
the constraints stability domain Dgt of (j4.ip ). The marginal in position variables of the distribution (j3.8p 
is the original distribution e~^^^'^'^dq which is independent of the mass-penalization v. 

Proof. The statement follows from reversibility and measure preserving properties of Verlet schemes 
(see [2]), from the Hybrid Monte Carlo rule ([HUM]), and from the construction of the mass-penalized 
Hamiltonian (Proposition I3.4p . □ 

4.3. Interpolation between exact and constrained dynamics. First we study the limit towards 
the leap frog/ Verlet scheme. 

Proposition 4.3 (SmaU penalty hmit). Assume that the IMMP integrator (|4.ip or (|4.2p is locally 
well-defined, and that u ^ 0. Then the latter two (computed with or without the Fixman correcting forces 
^qVfix^v) converge at order Oiv^) towards the leapfrog /Verlet scheme for the exact dynamics. 

Proof. Consider the auxiliary variables in the scaling {z,p^) — {^,^), and the following shift of 
Lagrange multipliers in (|4.2p . or equivalently in (|4.ip 

, 5t 

then the scheme (14.21) becomes 

P«+l/2 =Pn- y + %VHx,i/)(gn) + Vq^((7„) A„+i/2 

Pn+1/2 = Pn~ \i+l/2 
Qn+l ^ qn+ Vn+1/2 
Zn+1 = 2n + '^^A/7^Pj;+l/2 
£.{qn+l) = Zn+1 (C1/2) (4.4) 

Pn+1 = Pn+1/2 - y (VlC/ + V5Vfix,,.)(g„+l) + l^^yq^{qn+l)Xn+l 

Pn+1 = Pn+1 ~ "^n+l 

which converges at order 0{i'^) to a decoupled scheme where the {q,p) variables evolve according to the 
usual leapfrog/Verlet scheme, and the auxiliary variables {z,p^) are enforced by the constraints £_{q) — z. 
Now from (jl.3p . Vfi^.v is of order 0(i^^) which completes the proof. □ 

One can now construct schemes consistent with respect to the exact dynamics by letting the penalty v = 
vbt^ go to zero with the time step, for some fc > 0. Indeed, Proposition 13.61 shows that the mass-penalized 
dynamics (j3.5p converges towards the exact original dynamics for = at order 0(y'^\ Consequently most 
of the usual numerical schemes will be consistent at their own approximation order but bounded above by 
2k. The order of convergence refers to the maximal integer k such that the convergence of trajectories with 
respect to the uniform norm occurs at the rate 0{bt^^. Neglecting the order of the fluctuation/dissipation 
part (Step 2 of Scheme 14. ip , we deduce the following consistence property. 

Proposition 4.4 (Time-step consistency with exact dynamics). Assume that the integrator ()4.ip or 
(|4.2p are locally well-defined, and that v = VSt. Then the IMMP numerical scheme ()4.ip or (|4.2p (computed 
with or without the Fixman correcting forces 'S/qVf^^^^), is of the order 2 and consistent with respect to the 
original exact deterministic dynamics (j2.3p (i.e., with 7 = 0, a — Q). 
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Proof. Following the proof of Proposition 14.31 we have, by the implicit function theorem, locally, the 
RATTLE scheme is a standard leapfrog scheme (see [TH and the associated local mapping depends 
continuously on j/^ when v ^ 0. Therefore the usual calculation of the order of the leapfrog scheme (see 
[l4j ) holds uniformly with respect to v, and (|4.ip or (|4.2p are of order 2 consistent uniformly in v. Now, 
the IMMP Hamiltonian (|4.2p is a i^^ perturbation of the exact Hamiltonian, and the result then follows 
from applying a simple Gronwall argument. □ 

Similarly, we also obtain the limit towards the constrained/RATTLE scheme. 

Proposition 4.5 (Large penalty limit). Assume that the integrators (|4.ip or (14. 2p are locally well- 
defined, and that v — + +oo with ^ — > zq- Then the IMMP numerical scheme (14. ip or (j4.2p (computed with 
or without the Fixman correcting forces VgVfix.zyj converges at order -i^ towards the constrained/RATTLE 
scheme for the rigid constraints CC'Zn+i) = (also computed with or without the associated geometric 
correcting Fixman forces VgVfjx defined by (13. 3p ). 

Proof. Following similar steps as in the proof of Proposition 14.31 with auxiliary variables in the scaling 
{z,p^) = (^, 2-)^ and the following shift of Lagrange multiphers: 

we obtain convergence at the order to a decoupled scheme where the {q,p) variables evolve according 
to the usual constrained/RATTLE scheme. Now from (|3.3p . Vfix.i^ is a perturbation of Vnx of order 
which completes the proof. □ 

Remark 4.6. We conclude this section with a practical recipe for tuning the mass matrix penalty. 
Identifying a suitable value of v can be done, for instance, by computing the time averaged energy error 
(|5.ip in a dynamical simulation, or the Metropolis acceptance ratio (|5.2p in a sampling simulation, which 
gives a precise quantification of the time-step error. Then increasing the penalty v can save computational 
time as long as it leads to a reduction of the time-step error. Indeed, this means that the selected fDOFs 
are limiting the time-step stability region. Prescribing the time-step error, a maximal time-step Stmax 
associated with the largest penalty Vmax that is able to improve stability can be obtained in this way. 
Finally, one can set, for example, v = ^"""^ St in Scheme 14.11 to obtain an order two convergent scheme 
with an increased stability region. 

5. Numerical simulations for the A^-alkane model. The IMMP method is numerically tested on 
the united atom A^-alkane model. The united atom model is a coarse-grained description of linear alkane 
isomers in which the hydrogen atoms are not resolved and the molecule is modeled by a chain of particles 
interacting with effective potentials, see, e.g., [30]. The integrator studied in the present section is described 
in Section [3] and Section 31 and is systematically compared with the exact dynamics, which is numerically 
integrated by a simple Leapfrog/Verlet scheme with or without a thermostat, and with the constrained 
dynamics numerically integrated by the RATTLE scheme with or without a thermostat. Constraints are 
resolved using a simple Newton algorithm with a Gaussian linear solver. The Fixman forces are not resolved 
in the dynamical part, and thus the dynamics interpolates for a large penalty with the "rigidly" constrained 
dynamics. However, they are used in the Metropolis rule when sampling is considered. 

The Metropolis/HMC step for the dynamics integrator as described by Scheme l4l2l is added only when 
sampling is studied. 





Verlet = 0.5 j/ 1.0 = 1.3 RATTLE 




.024(1) .032 (1) .046 (1) .059(1) .077(1) .093 (2) 




.013(1) .014(1) .022 (1) .028(1) .035 (1) .049 (1) 



Table 5.1 

Critical time steps tc^" and 54^^™''' of the butane deterministic dynamics and of the sampling schemes, respectively. 
Note the interpolation property from constrained dynamics to exact dynamics. 
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Length Evo., Verlet 



Length Evo., nu=0.26 



Length Evo., nu=0.69 




1 23456789 10 

Physical Time 
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2.235- 
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Physical Time 



01 23456789 10 
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Physical Time 




1 23456789 10 
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2.245t 
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Figure 5.1. Oscillations of the butane end-to-end length, for the Verlet integrator, the IMMP integrator, and the 
constrained integrator. Simulations are computed with a prescribed initial energy. Note the interpolation property of the 
IMMP scheme. This figure is associated to the frequency analysis in Figure f5. 21 



5.1. The A^-alkane model. The model consists of a chain of iV atoms with position vectors g'*-* £ R^, 
hence the vector of DOFs is q — {q'^^\ . . . , q'^^^) G R'^^. The mass of particles is normalized to be mM^ = 1 
for all i = 1,...,N. The interaction potential V consists of three short-range potentials that involve 
2-body, 3-body, and 4-body terms 

V{q)= Vbo.M'\q^'^)+ E V^nsUQ^'\q^'\q^'^)+ ^ torsion (g« , -Z^^'^ 'Z*'^^ 9^'^ ) • 

In the presented simulations we focus on short-range interactions as they are primarily responsible for 
the stiffness of the system. Thus we omit the long-range Coulomb or Lcnard- Jones interaction potentials. 

2- body interactions. The pair-wise bond potential Vbond depends on the distance ~ \q^'^^^^ — g*-*-*], 
i = l,...,iV~lof the bonded particles (g*-*-* , c^'^^^'^) in the linear chain. Typically this potential is assumed 
to be harmonic with the equilibrium distance r^'-* = rg. In all simulations considered here the bonds are 
treated as rigid bonds of the constant length = 1. Thus defining the constrained DOFs 

aond(<?) = (rW,...,r(^-^)) = (l,...,l). 

3- hody interactions. The interaction of three consecutive particles (g^*-*, g*^*"*"^-*, g^*"'"^-') in the chain is defined 
by the potential 14ngic (fi*^*' ),i = l,...,iV— 2 that depends on the angle between the vectors (g^*+^^ — g(*+i) ) 
and (g'*+^^ — g'*^). In the simulated model we assume the bending angle potential 

Kngic(0) = ^sin2(0-|). 
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N=4, Frequency 
Cumulative Distribution 




Figure 5.2. Frequency distribution of the end-to-end length oscillations of the butane Verlet dynamics and of the 
IMMP dynamics, cf. Figure [5.11 Note the two main frequencies (slow torsion and fast bond angles oscillations), and the 
slow components due to fluttering of resonances. 

Depending on simulation the bending angles can be mass-penalized or directly constrained leading to the 
definition of penalized DOFs 

Cangio('i') = {()^^\ ■ ■ ■ , 9^^^'^^) —— in the case of mass-penalization, 
Cangicl?) = i9'"^\ ■ ■ ■ 1 0^^^^^) — in the case of rigid constraints, 
where Za G M^^^ are the auxiliary variables of the IMMP method. 

4-body interactions. The last contribution Vtorsion(0^*^) to the interaction potential depends on the torsion 
(dihedral) angles 0*^*), i = 1, ... — 3 which are defined as angles between two planes spanKq^^^'^) — 
(g(^+2) - and span{(q(*+2) - - q^'^)}. We use a simple choice of the torsion 

potential 

Viorsion(0) = -Bq COs{(j)) . 

The torsion DOFs are mass-penalized with the auxiliary variable zj, S R^"'^ 

System of units. For the purpose of computational tests we have chosen the system of units in which the 
mass of united atoms is m^^^ — ■ ■ ■ — m*^^-* — ijiq — 1, the equilibrium bond lengths are r^^-* = . • • = 
j-iN-i) — — 1 and the inverse temperature /3 = 1/kT = 1 at the ambient temperature T — 300 K. The 
time step 6t — 0.01 in these units corresponds to the physical time step 6t = ^/ /3moro = 3 is. Following 
physically relevant parameters (see, for example, [331 [50]) with an artificially slightly stiffer bending angle 
potential leads to the angle potential constant Aq = 500 and the torsion potential constant i?o = 20. 
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Figure 5.3. Interpolation of short trajectories of the butane length with IMMP dynamics, from constrained to Verlet 
dynamics. The marked trajectory corresponds to v = 5.5. The Verlet dynamics is the oscillating limit. 



Applying direct constraints to bonds or bending angles in the presented formulation of the IMMP 
method consists in replacing the auxiliary variables z associated with fully constrained DOFs by 0. As 
described in Proposition 13.71 on the large penalty limit (see also Appendix |A] and Section [6] for the limit 
of stiff distributions), the associated dynamics is the standard "rigid" constrained dynamics. Then for 
sampling, the Fixman geometric corrector (II. 4|) is computed and accounts for the geometric bias of all the 
constrained DOFs including bonds. Note that the physical relevance of the geometric bias due to bonds 
is ambiguous since the quantum resolution may need to be introduced for the bond interactions. Such a 
change can be easily incorporated by computing the Fixman corrector potential associated with the bonds 
only and subtracting it from the Fixman corrector potential related to all the constrained DOFs. 

5.2. Efficiency criteria. The efficiency of the IMMP method as compared to the Verlet /Leapfrog 
integrator is quantified from two different viewpoints. 

Dynamics. To compare computational efficiency of numerical schemes, a notion of critical time step has 
to be introduced. The critical time step JtjJ^" is defined implicitly through the formula 

[H^+l-H^]+^l,{pZ,q^)=a, (5.1) 

where Nc denotes the number of rigidly constrained DOFs, [•]+ is the positive part, is the energy of the 
numerical scheme at the n-th step, and /i^ is the Boltzmann distribution associated with the Hamiltonian 
at hand as defined in l|3.8p . The quantity a is a prescribed typical non-dimensional error of the energy 
per degrees of freedom as compared to the inverse temperature (3. At least when a is small, this defines 
uniquely (5tjJ^". Thus for a given a, the larger J^jJ^" is, the less costly the method is in terms of force 
evaluations per integration time step. 

The critical time step Sf^^^ is compared for the IMMP method and the Verlet/Leapfrog scheme. To 
achieve a fair comparison, it has then to be checked that the IMMP penalty, which introduces a dynamical 
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Figure 5.4. Convergence in terms of the pathwise error of the short time butane length IMMP dynamics compared 
to Verlet dynamics and constrained dynamics, with orders of convergence 0{u'^) and 0{u^^). 



modification, does not modify the relevant slow frequencies and slow varying components of the system 
on large time intervals. This can be done by comparing the distortion of the frequency distribution of a 
long trajectory. In the numerical tests we analyze for t G [0, T] the trajectory of the end-to-end molecule 
length t I—!- L{t), and we define the frequency density as the normalized square modulus of its Fourier 
transform d{uj) = \L(uj)\'^ /Z f where Zf is the normalization making d(uj) a probability density function, 
and L{lu) = L{t) exp{iujt) dt. We plot the cumulative distribution function 



D{lo) 



d{uj')duj' . 



Hence, in the figures dominant frequencies correspond to jumps in the associated cumulative distribution. 
Of course, this comparison remains still largely qualitative. 

Sampling. When using a Metropolis correction step in the numerical simulation, as in Scheme 14.21 the 
critical time step is correlated to the rejection rate of the Metropolis step, since the larger the former is, 
the more rejections will occur. Thus for sampling methods, the time step 5t^^'"P' is tuned in order to 
achieve a given rate of rejection p satisfying 



j exp (-/3[i/„+i - i/„]+) g") = 1 - P, 



(5.2) 



where exp (— — Hn\^) is the Metropolis weight that appears in Generalized Hybrid Monte-Carlo 
methods, as defined in Scheme 14.21 Then, computational efficiency is defined by comparing (5t^^™P' with 
the mixing time in terms of physical time of the resulting Markov Chain. Equivalently, and more directly, 
the results are presented as the mixing time in terms of iteration steps. Such comparison requires choosing 
an appropriate notion of the mixing time. At least when the momenta are overdamped (i.e., 7 — > -l-oo), 
the chain becomes reversible, and the mixing time can be rigorously defined using the spectral gap of the 
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Figure 5.5. Left: Equilibrium PDF of the end-to-end length of the butane molecule with the GHMC scheme, using 
Verlet, IMMP (penalty v), and constrained/RATTLE integrators. Note that the constrained integrator does not sample the 
correct measure. Right: The autocorrelation function in terms of iteration steps for the GHMC scheme comparing the IMMP 
and the Verlet integrator. The decrease in the -decorrelation time 15.3(1 is by the factor 1.8. 



underlying Markov kernel. In the present work, the decorrelation time of relevant observables is used. 
Assuming the initial state qo of the system is at equilibrium, the normalized autocorrelation function 
associated with a given position observable q i— > is 

_ E[0(g„)0(go)] 
^"^^^ - mqof] ' 

and it can be computed in large time simulations using path averages. Then ^^-decorrelation time ncorr, 
calculated in terms of iteration steps, is given by the formula 

+ OC 



2^C„(0)2. (5.3) 



n=0 

This quantity corresponds to the approximate number of steps of the chain for the observable 4> to decor- 
relate from its past values. 

5.3. Numerical results. For comparisons we choose an often studied observable defined as the end- 
to-end distance of the alkane chain. 

I. The butane model {N = 4). The bonds between atoms are rigidly constrained, and bond angles are 
treated using the IMMP method with Scheme \4A[ 

Dynamical behavior. The dynamics of the butane model is studied using deterministic dynamics with a 
prescribed initial energy. As shown in Figure EHJ the IMMP dynamics with different penalty i/ yield an 
interpolation from the mixed torsion/bond angles oscillations of the exact dynamics, to the simple torsion 
oscillation of the constrained dynamics. Depending on the frequency introduced by the mass penalization, 
some fluttering resonance can be observed in Figure 15.11 between the torsion and the bond angles. 
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Behavior observed in Figure [?TT] is related to the frequency density of the end-to-end length oscillations 
which are depicted in Figure 15.21 The bond angles oscillation frequency appears to slow down with the 
IMMP penalization and is eventually setting to a single frequency of the constrained dynamics. Note the 
slow modes introduced by the resonances. 

The dynamical interpolation of the IMMP from the exact dynamics to the constrained dynamics is 
demonstrated on short time trajectories in Figure [5T51 The associated convergence orders, 0{v^) and 
0{i'~'^), respectively, are captured in Figure WM The time step stability is studied in Table [5T] The 
IMMP method enables an increase of the critical time step of the Verlet scheme, following the interpolation 
property. 

Sampling behavior. Exact sampling of the equilibrium distribution on a very large time scale, whatever 
the value of the IMMP penalization, is shown in Figure 15.51 The distribution of the butane length for 
constrained bond angles is clearly distorted. The mixing time to equilibrium is also studied. The autocor- 
relation function of the length evolution in terms of iteration steps is plotted in Figure [531 the faster 
convergence of the IMMP method is demonstrated. The £^-decorrelation time (|5.3p is decreased by the 
factor 1.8 using the IMMP method. 




Figure 5.6. The trajectory of the N = 15-alkane dynamics for the Verlet scheme and the IMMP scheme. Note that 
the IMMP penalty does not modify substantially the low frequencies/slowly varying components. The frequency analysis is 
presented in Figure \57f\ 



II. The ALKANE model {N = 5, . . . , 20). In this test the bonds and bending angles between atoms are 
rigidly constrained, and torsion angles are treated using the IMMP method. Note that the rigid constraints 
on torsions would lead to a rigid molecule, losing completely the evolution of the molecule length. The 
IMMP penalty is increased with the system size using the linear scaling v = vN . In Section[7]we present a 
mathematical justification of this scaling for linear systems. In principle, introducing inertia in the torsion 
angles adds weights on the diagonal of the mass-tensor in the internal coordinates of the molecule. This can 
reduce the lowest eigenvalues of the mass-tensor, and can explain the reduction of the multi-scale nature 
of the system that are observed in the simulations below. 
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Figure 5.7. Spectral densities of the end-to-end length of the alkane for N ■ 
does not modify substantially the low frequencies/slowly varying components. 



■ 5, 10, 15, 20. Note that the IMMP penalty 



Dynamical behavior. The frequency of oscillations of the alkane length in Figure \5M and in the spectral 
analysis in Figure 15.71 are not substantially modified by the IMMP penalization. One can observe a small 
group of fast oscillations in the middle of the Verlet dynamics plot in Figure 15.61 which is not present in 
the IMMP case. This translates in the top of the spectral plots in Figure 15.71 where a cut-off of the fastest 
oscillatory scales for the IMMP case occurs. 

The critical time steps (5tJ?^" with respect to the system size are depicted in Figure [5781 The gain in 
time stepping increases more than linearly with the system size N . This behavior, however, depends on 
the initial scaling i' = vN of the IMMP penalty. 

Sampling behavior. Exact sampling of the equilibrium distribution on very large time scales (whatever the 
value of the IMMP penalization) is shown for TV = 10 in Figure [5T9l with the auto-correlation of the latter 
observable with the gain in mixing time of the IMMP dynamics. The precise ratio of the £^-decorrelation 
time between the IMMP integrator and the Verlet one is given in Figure [5. Ill for different system sizes. It 
increases again more than linearly in iV, in fact exponentially for this particular system. The associated 
critical time steps i5t^'^™P' are in Figure [5.101 We observe that the critical time step increases with large 
N which demonstrates that in the present case, the IMMP method heals the decrease of the Metropolis 
rejection rate of for large systems (see also |24) ) . 

Conclusions. The presented numerical studies demonstrate that for integrating the dynamics, the IMMP 
allows for relaxation of the time-step stability restrictions. In the case of sampling methods the IMMP 
method decreases the decorrelation time measured in terms of Monte-Carlo iteration steps leading to more 
efficient sampling algorithm. In both cases the improvements are increasing with system size N . 

6. The infinite stifi"ness limit. Throughout this section, one introduces a potential function with an 
explicit dependence with respect to the fast variables (9,2;) ^ U{q,z) together with a stiffness parameter 
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Dynamics: Critical 
Time Step 
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Figure 5.8. Critical time steps, with error bars, (Stc^" of the Verlet and IMMP dynamics with respect to the system size. 



e. The potential energy V can then be written in the form 

V{q) = U{q, ) , with a confining assumption inf U{q, z) > K{z) , 



where K{z) > clog|z| as \z\ oo. The functions ^ are then indeed "fast" degrees of freedom (fDOFs) in 
the Umit e ^ 0, the system being confined to the slow sub-manifold Mo = {? I ^(9) = 0}. 

We prove, that under appropriate scaling of the mass penalty i^^ = the IMMP method is asymp- 
totically stable in the stiff limit, converging towards standard effective dynamics on the slow manifold 
Mo. 



6.1. Thermostatted stiff systems. The canonical distribution becomes 



Z. 



(6.1) 



In the infinite stiffness limit (e 0) the measure concentrates on the slow manifold A^o- The limit is 
computed using the co-area formula (see Appendix [Xl for relevant definitions of surface measures). In order 
to characterize the limiting measure we introduce the effective potential 



Kff(<?)--^ln J e-^^(''^) dz. 



(6.2) 



Lemma 6.1. In the infinite stiffness limit (e Q), the highly oscillatory canonical distribution (j6.ip 
converges /Lte^/xo (in distribution) towards ^,o{dpdq), which is supported on Mo, and defined as 



Mdpdq) = ^e-/3(^p"M-V+^4«(<?)) dpS^^,)=o{dq) 
^0 



(6.3) 
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N=10 Alkane Length N=10 Autocorrelation 

Equilibrium PDF Function 




Iteration Steps 

Figure 5.9. Left: Equilibrium PDF of the end-to-end length of the alkane chain with the GHMC scheme, using Verlet 
and IMMP integrator (penalty v). Right: The autocorrelation function in terms of iteration steps for the GHMC scheme 
comparing the IMMP and the Verlet integrator. The system size is N = 10 

Its marginal distribution in position is given, up to the normalization, by 

e-PyM,)s^^^)^„idq). (6.4) 

Proof. It is sufficient to consider distributions in the position variable q only. Let be a S- 
neighborhood of Aig where dq = ^'^5^(q)=tz{dq) dz. We construct a decomposition ip = 'pi-\-ip2 of continuous 
bounded observables such that supp(^i C and supp(/52 n W'/^ = 0. Using the confining property of 
U{q,-) we obtain 

By continuity of e i-^ J V'i('z)6~'"^''*'^''<^4(g)=ez('^9) ^^id by the dominated convergence theorem 

and the result follows after normalization. □ 

The infinite stiffness limit (e 0) of highly oscillatory dynamics has been studied in a series of papers 
[37l [39l [27l [5l [341 [35] . The limiting dynamics can be fully characterized in special cases. For example, when 
the highly oscillatory potential is linear and non-resonant (at least almost everywhere on the trajectory, 
see [39]), it can be described through adiabatic effective potentials. See also [lOl [6] for a recent work on 
some related numerical issues. However, when the system is thermostatted, one can postulate an "ad hoc" 
effective dynamics (|35|) exhibiting the appropriate limiting canonical distribution given by (j6.4p . Such 
dynamics can be obtained by constraining the system to the slow manifold A4o, and adding a correcting 
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Sampling: Critical Time Step 
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nu =0.12(N-3) 
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Figure 5.10. Critical time steps, with error bars, St^'^'^'^ of the Verlet and IMMP dynamics with respect to the system 
size. Note that the IMMP penalty heals the degeneracy of the rejection rate for large systems. 



entropic potential, sometimes called Fixman corrector from [T^, which is due to the geometry of Aloi ^^nd 
is given by 

14ix(g) = ^ln(detG(g)) , (6.5) 

where G{q) is the n x n Gram matrix defined in p.2p . 

In general, since the effective potential (|6.2p is not explicit, one may need to couple the system with 
virtual fast degrees of freedom to enforce the appropriate effective dynamics associated with (|6.2p . The 
resulting extended Hamiltonian is then defined on the state space T* {Mo x K") (the cotangent bundle) 
and is given by 



ifeff (p,P., g, z) = ^p'^M-^p + ^pIm^^p, + U{q, z) + yfi,(g) 

e('z)-o. 



(6.6) 



Definition 6.2 (Effective Langevin process with constraints). The constrained Langevin process 
associated with Hamiltonian (16. 6|) is defined by the following stochastic differential equations 



(6.7) 



' q = M^V 
z = M-V. 

p ^ -ViUiq, z) - VqVa^iq) --fq + aW- Va^ A 
Pz = -V2C/ (g, z) - jzz + a^Wz 
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L2 Decorrelation Time 
Leapfrog/Verlet vs IMMP Ratio 
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N = System Size 

Figure 5.11. Comparison the -decorrelation time of the end-to-end alkane length in terms of Monte- Carlo iteration 
steps. The ratio between the case of the Verlet integration and the case of the IMMP integration is depicted. Note that the 
vertical axis is in the logarithmic scale and an exponential gain occurs. Error bars are plotted, but do not appear at this 
scale. 

where Vi and V2 are respectively derivatives with respect to the first and second variable of the function 
U{q,z), W (resp. Wz ) is the standard multi- dimensional white noise, 7 (resp. ^z) a dx d (resp. n x n) 
symmetric positive semi-definite dissipation matrix, a (resp. Cz) is the fluctuation matrix satisfying aa^ = 
^7 (resp. Oz<^'^ = ^^^)- -^^^ processes A G M" are Lagrange multipliers associated with the constraints 
(C) and adapted with respect to the white noise. 

We formulate reversibility of this process as a separate lemma. 

Lemma 6.3. The process defined in ^6. 7p is reversible with respect to the associated canonical distri- 
bution whose marginal distribution in {q,p) variables is 

^eff 

with the q-marginal 

When 7 and 7^ are strictly positive definite, the process is ergodic. 

Proof. The process (|6.7p is a Langevin process with mechanical constraints, exhibiting reversibility 
properties with respect to the associated Boltzmann canonical measure (see the summary in Appendix |B|) . 
Then the q-marginal is obtained by remarking that the integration of any function of ^p"'" M~^p-\-^p^ M~^pz 
with respect to dpz cft* Moi.dp) results in a constant independent of q. □ 

The properties of thermostatted highly oscillatory systems are summarized in Table WT\ 
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Finite stiffness 


Infinite stiffness limit 


Infinite stiffness 




e > 





e = 


Dynamics 


Highly 


Adiabatic 


Effective with 




oscillatory 


(if non-resonant) 


constraints 




+ fluct./diss. 


+ non-Markov fluct./diss. 


+ fluct./diss. 


Statistics 


Canonical 


Positions on A^O: 
free velocities. 


Canonical on T*Mo, 
geometric corrector. 


Numerics 


Leapfrog / Verlet 


Time-step 


Lcapfrog/Verlet with 




+ fluct./diss. 


restrictions {6t = o(e)) 


constraints 
-1- fluct./diss. 



Table 6.1 

Stiff Hamiltonian systems and associated commonly used numerical methods (Mo denotes the slow manifold). Two 
different schemes are required for the stiff system and its effective Markovian approximation. 



6.2. Stability of the IMMP dynamics. We assume that the mass-matrix penalty parameter i/ = i/g 
grows to infinity in such a way that linie^o ev^^ = v. 

The original Hamiltonian with the stiffness parameter is expressed explicitly as 

7J,(p,g) = lp^M-V + [/(g,^), (6.9) 
and including the mass-matrix penalization one gets 

H.APu.,q) = \pIM-^P.^ + U{q, ^) -f Vii^,,M , (6.10) 
or in its implicit formulation 

11 z 
HiMMp{q,z,p,p^) = -p'^M'^p+ -p^M^^p^_ + U{q, — ) + Vfix,i.,(g) , 
2 2 j/ge 

aq)^-z- (a J 



(6.11) 



One immediately sees that i?iMMP is non-singular when e — > and converges to the effective Hamiltonian 
on the slow manifold, 

H,sAq^z,p,p,) = ip3^M~V+ ^P^^CV + U{q, |) + VM ^2) 
e(g) = 0. (C). 

The expression (|6.11[) represents a minor generalization of ffcff in (|6.6p , but it leads to the same canonical 
marginal distribution fLas{dpdq) in (p, q) variables as given by (|6.8|1 . The continuity in e of Himmv implies 
stability of the associated dynamics and their numerical integrators. We first derive the limits of the 
original and penalized canonical distribution. 

Proposition 6.4 (Limits of canonical distributions). Consider the canonical distributions fj,i,^{dpdq 
associated with the mass penalized Hamiltonian (|6.10[) . but considered with respect to the variables (p = 
MM~^p^^,q). In the sense of weak convergence of measures we have fii,^^ficS as e Q with /ieS defined 
by (Ell). ' 

Proof. The first convergence is proved in Lemma l6. II For the second one, the following notation will 
be used 

SqMdq) = S^{q)=czidq) , and Sp^^p^{dp) = S^t M-ivi{q)=eMr'pMP'> ■ 
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To prove the convergence towards /lofr we consider a (5-neighborhood lA^ of where 

dpdq = . ^ 6q,ezidq) dz dp,ep, (dp) dp, , 
det M2 

and a decomposition of the bounded observable (in (p, q) variables) ip = + such that supp fi C 
and supp(y92 H W'/^ = 0. Thus, keeping in mind that p^^ = M^^M~^p, and using the confining property of 
the potential U{q, •) we obtain 

J >p{p,q)e-^""^dp,^ dq = J ipi{p,q)e-f^""'dp,^ + ©(e'^^^^/^)) = + ©(g-'^^t^/^^)) . (6.13) 

Applying the change of variables p,^^ — Mi,^M~-^p yields 

dp^ = det (Mj, Af-i) dp = 1/2" det det (G + \m~^) dp , 

and setting eM~^pz = V^^ M^^p and = ^(g) we get 

1 z 

Hu,{Pv,,q) = T^P^ M~'^p + vle^pl M~'^p, + t/(q, ) + Va^.^M) = HiuMpiq, z,p,Pz) ■ 

2 i^^e 

Thus substituting back to (|6.13p we obtain 

= (i/ee)'"y" V'le-'^^^^^^^'^'P'P^Met {G + ^M-^) Sp^.pjdp) dp, Sg^dq) dz , 

and thus 

/, — ^ y ^^e-''^''«.-(«'^'f'f^)det (G) (dp) dp, <5^(,)=o(rf9) dz . 
Using the co-area formula we obtain 

det (G) ^vc(g)J\/-ip=o(c?P)'5e(g)=o(c?'7) = <^T'M„{dpdq) , 

which leads to the final result after integration of the (p,, z) variables and normalization. □ 

Remark 6.5. Due to the fast oscillations, the distribution of impulses in the limiting distribution 
fxo in (|6.3p is uncorrelated, whereas after the mass-matrix penalization, the limiting distribution (|3.8p has 
almost surely co-tangent impulses (i.e., satisfying the constraints Vq^M~^p = 0). This explains the role 
of the corrected potential energy Vbx taking into account the curvature of A^o- 

In the next step we inspect the infinite stiffness asymptotic of the penalized dynamics. 
Proposition 6.6 (Infinite stiflFness limit). When e ^ with v = v^ ^ ^ and V{q,^{q)) = U{q, 
the IMMP Langevin stochastic process (13. 5p converges weakly towards the following coupled limiting pro- 
cesses with constraints 

'q^M-'p, 
p = ~ViU{q, |) - VqFfix(g) --iq + aW ~ , 

C(g) = 0, (G) (6.14) 

i = M^^Pz , 

1 z 
Pz = --V2C/(g, -) - -fzZ + (JzW, . 
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where Vi and V2 are respectively derivatives with respect to the first and second variable of the function 
U{q,z), and {Xt\t>Q adapted stochastic processes defining the Lagrange multipliers associated with the 
constraints (C). 

The process {qt,Pt\t>o defines an effective dynamics with constraints fDefinition \6. 2\) for thermostatted 
highly oscillatory systems. It is reversible with respect to its stationary canonical distribution given by /ioff 
(|6.8[) . and is ergodic when (7,7z) are strictly positive definite. 

Proof. The proof is similar to the proof of Proposition 13.71 Here we have 

and (|3.5p translates, up to a change of Lagrange multipliers, into 

' q = M-^p 
z ^ M-^p, 

P = -ViC/ - VqVfi^.^(g) - 79 + crM^ - V,^ A 

1 . A ^^-^^^ 

Pz = V2[/ - + a-zWz H 

= (C.J. 

The rest follows the proof of Proposition [3171 
□ 

Remark 6.7. When v +00, by a classical averaging argument (see, e.g., [27]), one can check that 
the limiting dynamics are the effective dynamics pointed out in [35] 

q = M-^p 

P = -VgUMq) - V,14x ~-/q + c7W- A (6.16) 

with the stationary canonical distribution (|6.8p . 

6.3. Stability of the IMMP integrator. The numerical scheme (Scheme 14.11 proposed for the 
IMMP method p.Sp ) is also stable in the limit of infinite stiffness e ^ 0. Recall that we consider a 
reversible, measure preserving numerical flow ^gl{p,Pz,q, z) associated with Hamiltonian (|6.1ip Himmp 
with constraints (modified potentials could similarly be considered). 

Proposition 6.8 (Asymptotic stability). In the limit ev^ v, the numerical flow (f>^j associated 
with the leapfrog /Verlet integrator with constraints for the IMMP Hamiltonian (|6.1ip converges towards 
the numerical flow <i>^(, which is the leapfrog /Verlet integrator with geometric constraints associated with 
effective Hamiltonian (|6.12p on the slow manifold. 

Proof. The statement is a direct consequence of the implicit function theorem and the continuity of 
the leapfrog integrator with constraints (|4.ip with respect to the parameter P = eVe- Indeed, considering 
the shift of Lagrange multipliers A A + iV2J7 and taking the limit e ^ we obtain the appropriate 
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leapfrog scheme 



Pn+l/2 Pn - y ViC/(g„, ^) - Vq^(g„)A„+i/2 



5t 



v 



Pn+1/2 =Pn - ^V2f/(q„, 
= gn + (5iAf~Vri+l/2 

= 

Pn+1 = Pn+1/2 - yVi?7((7„+i, f!^) _ Vq^(g„+i)A„+i 

2 Z '^^ 7-7/ Zn+1 X 

=P«+1 - Tp^2U{qn+l,^^) 

I V,e(g«+i)M"Vn+i = 0. 



(C1/2) 



{Gi) 



□ 

By convergence of the Hamiltonian (|6.1ip to (|6.12p , similar asymptotic stability properties holds when 
a Metropohs step is introduced. 

The results and properties discussed in this section are summarized in Table 16.21 





Zero mass 


Positive 


Infinite 




penalization 


mass-penalization 


stiffness limit 






e, > 




Dynamics 


Highly oscillatory 


IMMP 


Effective with 




+ fluct./diss. 


+ fluct./diss. 


constraints+ fluct./diss. 


Statistics 


Canonical 


Canonical with 
correlated velocities 


Canonical on 


Numerics 


IMMP + fluct./diss. 



Table 6.2 



The IMMP dynamics and the Verlet numerical integration are both asymptotically stable in the infinite stiffness regime 
if ^ P < +00. If the mass-penalization vanishes (v = 0) one recovers the original physical stiff system. The canonical 
distribution is always exact in the position variable. Notice that due to the penalized mass-matrix (v > 0) the statistics have 
correlated velocities. 



7. Numerical analysis of a harmonic particle chain. In this section we present rigorous analysis 
for a special case of the linear chain with harmonic interactions. The analysis supports scaling properties 
of the IMMP, with respect to the size of the chain, observed in numerical simulations of the general linear 
alkane chains. We consider the thermodynamic limit N +00 where TV is the size of the system. It is 
shown that the macroscopic dynamics of the IMMP method behaves continuously (uniformly with N, and 
in the L2 norm for the position profile) with respect to the re-scaled mass penalty parameter 9. 

At the same time, the time-step stability of the IMMP numerical scheme (|4.3I) is compared with the 
standard Verlet scheme, and the critical time step is shown to be increased by a factor 9N. 

From the spectral point of view, the IMMP method behaves in this linear case as a low-pass filter. This 
proves, in this simplified case, the ability of IMMP method to respect macroscopic dynamical equivalence, 
while saving computational time up to a factor of order 0{N). 

7.1. Conservation of macroscopic dynamics. The model we consider consists of a chain of par- 
ticles which interact through the harmonic (quadratic) potential Wint('') = r'^/'2. Each particle is also 
individually submitted to a macroscopic confining exterior (quadratic) potential Woxt('")- After converting 
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to the non-dimensional form the typical quantities involved in the model enable us to write a scaling at 
the mass-transport level where the dynamics of the chain is described by the Hamiltonian 



1 



T 
P P- 



N-1 

E 

1=1 



N 

E 

1=1 



(7.1) 



and by a coupling with an exterior thermal bath at the re-scaled inverse temperature (3^ = I3N~^. In the 
expression (|7.ip the functions r g M Wint {t) S M and q ^ Wext (?) S K are the smooth interaction potential 
and the exterior potential, respectively. The linear operator : M^^^, having the components 



qt+i 



l/N 



1, 



,7V- 1. 



represents the discrete gradient associated to the chain with the Neumann boundary conditions. Its 
transpose operator is denoted (V*)"^ : M^"^ K.^. The discrete Laplace operator is then defined as 
Arf = — (V'')"^V''. The particles are represented by their re-scaled positions q = {qi, ...,qN), so that the 
typical position and deviation of q is formally of order 1 with respect to TV. This can be seen by considering 
particles in the chain as indexed by a; = G [0, 1]. We choose to work with such scaling in N that it 
prescribes the macroscopic timescale of the chain profile at order one with respect to N. 
Following our general construction we obtain the mass-penalized Hamiltonian 



1 



pI^^^i^^Pi^n 



N 

E 

1=1 



N 



^'int(Vfg) -f ^i;cxt(9») 



(7.2) 



i=l 



We chose the penalizing matrix to be the identity matrix A/^ — Id , hence the penalized mass-tensor 
becomes M^^^ = Id — v'^ l^d, and the fluctuation/dissipation tensor is taken proportional to the identity 
matrix. The system of stochastically perturbed equations of motions then becomes 



■ <7 = (Id - u'^^dy^p 
p = Arf (7 - (g) - 79 + (7 VA^IF , 



(7.3) 



with fluctuation/dissipation identity cr^ = 2j3~^^. The associated canonical equilibrium distribution is 
then given by the re-scaled inverse temperature (3n = (3N~^. 

In order to treat the limit iV — > 00, we introduce the ^2-norm in the position space 



as well as the /i_i-norm in the momentum space 



1 ^ 1 

1=1 



IbllL, 



N 2 . 

2=1 £^ \ i — 1 



In the above expression, jj- X]i=i Pi <^an be seen as the orthogonal projection in £^ on the one dimensional 
kernel of the Neumann discrete Laplacian A^. The quadratic form ||g||Q + endows the phase-space 

with a Hilbert space structure. 

Proposition 7.1 (Convergence of the macroscopic dynamics). Assume that the exterior potential 
Vcxt is bounded and that its derivative satisfies the Lipschitz condition 



\WcKtiQ2 



(gi)lk_, < Ly \\q2 - qiWi , 
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where Ly is independent of N. For any T > 0, let t {p'^ (t) , q'^ (t)) be the solution, for t e [0,T], of the 
evolution equation ^7.8^ with the initial condition 

(/(0),g'^(0)) = (AC;/V(0),g"(0)), 

where (p'^ (0) , q'^ (0)) is distributed according to the original equilibrium canonical distribution (associated 
with (|7.ip and /3jv — (3N^^). Then for all t G [0,T] one has the uniform convergence 



lim lim sup E 



\q'^{t)-q'^=^(t)\ 



Proof. We write X = {q,p), and introduce the norm 



0. 



\X\ 



The system (|7.3p becomes a stochastic differential equation in the form 

dXl = A^X^ + F{X^) + Y^dWt , 

where by definition 



(7.4) 



(Id - p2A^)-i 
Aw 



F{X) 









Na 



Duhamel formula gives an implicit expression for differences of solutions of (17. 4p with the same noise 



xr - X. 



iF{x^)ds + j:dWs) 



+e'^'{X^-X",)+ / e^^<^'-^\F{X:)~F{X'^))ds. 



(7.5) 



We estimate the individual terms on the right hand side in (|7.5p . We define P as the coordinate transfor- 
mation associated with the orthonormal spectral decomposition 



^Ad = p-Miag(5o 



, . . . , ojv-i 



where PP^ — Id . The eigenvalues of the discrete Neumann Laplacian are given, for A: = 0, . . . , iV — 1, by 



5*. = sin^ — 



2N J JV-*oo 



I 2 2 

k n 



(7.6) 



Denoting the spectral coordinates 



we have 



X^{q,p)^iN-'/^Pp,N-'^-'Pq) 



-1/2, 



N-1 



N-1 



k=l " fe=0 



The spectral decomposition leads to a block diagonal form of the operator e^"* with diagonal 2x2 blocks 
in the spectral basis 

ll ' 
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as well as for fc = 1, . . . , — 1 



6k 







where Ajy is the 2x2 block associated with the coordinates [qk,Pk)- Since Ac conserves the k-mode 
energy 8k{qkY + (1 + f'^<5fe)~^(pfc)^, one can check that for any > 1 the operator norm 

III e^"* mi < 2 + 2^2. 

— 1/2 

Similarly, since in the sense of symmetric matrices M^j^ < Id , we have the bound 

\\F{X'^)~F{Xm< M-^/^vWl')-<A<f)+lip'-p')) ^ <{Lp+^)\\X'^XX. 



h-i 



Using independence of Brownian increments we compute 



E 



Jo ^ v\ Jo II 



ds . 



Applying Gronwall lemma in (|7.5p and collecting all terms we obtain 



E 



\xr ~X?f_ 



< 



\x^,^xnl + ^T 



(7.7) 



where Ct is independent of N , and with X'^ being distributed canonically is given 



mT = sup 
te[o,T] 



N 



For a given random vector X such that 
imply 



+ E 



1X1 



< +0O, Parscval identity and the inequahty ||-||p < ||-||g 



E 



N-l 



«*)X||J] = ^E ||(e^-*-e^«*)l|[_ <2^ 



fc=i 



k=l 



X 



k,0 



(7.8) 



where ||-||j, p is the restriction to the k-th mode {qk,Pk)- Then one has, by orthogonality of P, 



i=l 







2 













1 ^ cr^ 

' 5k ~ 6k 



Up to normalization, the distribution of X° has the density e~ JV 2I]i=i "^xtCgi ) .^^j^^j^ aspect to the Gaussian 
distribution with the covariance matrix /3^^Id for momenta variables, and the covariance matrix (/JA^)"^ 
for positions. Thus we have the bound 



as well as 



lim 

N^oc 



E 



F{XO) 







2 




XO 








k,0 



< 2e''^l''''=="l'' 



SkP 



k.O 



< e 



4/3||«ext|| 



= E 



k.O 
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where T denotes the Fourier series expansion on [0, 1] with Neumann conditions, and (G*'fc)A;>i are canonical 

1 



centered Gaussian i.i.d. variables with the covariancc matrix /3 ^ 
the series is bounded 







By the Lipschitz assumption 



+ 00 



Since limi,^o \& * ~ e'^" *||| = 0, one can take the limit iV +oo and use the uniform convergence of 
the series in ()7.8|) to obtain limp^o liniA'^+oo ™t = in (|7.7p . The convergence of the initial condition 
limp^o limAT ^\Xl-Xl\_ follows by using similar arguments. The proof is complete. □ 

7.2. Relaxation of time-step stability restriction. To demonstrate improved stability properties 
of time integration algorithms we consider the IMMP scheme (|4.1|1 associated with the mass-matrix penal- 
ized Hamiltonian (|7.2p . Note that when the constraints are linear, the leapfrog scheme (RATTLE) applied 
to an implicit Hamiltonian is identical to the usual leapfrog scheme for the associated explicit Hamiltonian 
(|7.2p . We restrict the rigorous analysis to the quadratic interaction potential (fint(^) — ^^ro exterior 
potential (wcxt = 0); f^nd to the mass-matrix penalization operator (Id — The leapfrog scheme is 

defined as 

Pn+l/2 = + y (-Ad)g„ 

Pn+1 = Pn+1/2 + iry-'^djqn+l ■ 



k,0 



< (L, +7)E 







2" 




GO 








0_ 



fe=i 



. - (fc) 

=^0 t|| 



Denoting the spectral variables for k — I, 



,7V- 1 





(r 


4 \ 














Sk J 



1/2 



1/2 



NPp 
NPq, 



(7.9) 



we write 



where 



Lh — 



1 



n+l 1 _ r, 

n+l 



-hk + 
1 



hi' 



2 



and hk — St- 



1/2 



Since det (Lfe) — 1, the standard CFL stability condition is equivalent to 

\TT{Lk)\ <2 

which is fulfilled if and only if /ifc < 2 for all A: < — 1. Thus we arrive at the following bound on the time 

step 



St < 2 min 

0<fc<A' 



Sk 



1/2 
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Summarizing the above calculations and recalling (|7.6p we have the following characterization of the 
stability properties. 

Proposition 7.2. Suppose Vcxt — and consider a harmonic interaction potential Vint{r) — ^ with 
the mass-matrix penalization M^^ = Id — i>^lS.d. The leapfrog /Verlet integration of the IMMP harmonic 
Hamiltonian (|7.ip is stable in the spectral sense if and only if 



( \ 



1/2 



5t < 



sin 



V \ 2N J J 



{N ~ l)7r 



(7.10) 



Since we work with a Metropolis correction of the hybrid Monte-Carlo type, we are also interested in 
the limiting behavior of the energy variation compared to the temperature, i.e., 

f3N{H{pn+l,qn+l) - H{pn,qn)) , 

when (pmQn) are distributed according to the canonical distribution. This quantity gives the average 
acceptance rate of the Metropolis correction. The result we present here is similar to [4] where the authors 
analyze infinite dimensional sampling with the standard Metropolis-Hastings Markov chains. 

Proposition 7.3. Suppose Vg^t = o.i^'d consider a harmonic interaction potential Vintir) = ^ with the 
mass-matrix penalization M^,^ = Id — v^Ad- Suppose the state variable X = {Pun-^ i) ^■^ random variable 
distributed according to the canonical distribution associated with the mass-matrix penalized Hamiltonian 
()7.ip . Then the energy variation (3nAH after one step of the leapfrog integration scheme converges in 
distribution, up to normalization and centering, to the Gaussian random variable 

Law . /-/^ ^ 

>yv(0, 1) , 

with the mean and variance in the infinite size asymptotics for the IMMP method P > and St = St^ = o(l) 

N6t% , 2 N6t% 
n^N ~ ^ , ana Gat ^ ^ , 

and for the Verlet integration of exact dynamics with 6t = dt^ = o{l/N) 

5 5 
~, o^'^^^N > and a% - T^^'^'^^ ■ 

Proof. We start with a canonically distributed state X — {q,p), which is, by assumption on the form 
of the interaction potential, a Gaussian random vector. After changing to the spectral coordinates (|7.9p 
we have the spectral representation of the Hamiltonian 

N-l rl/2 
k=l ^ 

and introducing Gaussian random vectors IJ and V with the identity covariance matrix we can write 
We then compute explicitly the change of the Hamiltonian after one step of the leapfrog integration 



Pr^AH = J2 i (^^ ^ {LlL, - Id ) (j;^- ) . (7.11) 
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Since det (L'^L^) — 1 the matrix Lj^L^ — Id has two positive eigenvalues (A^, — 1, l/X^ — 1) which satisfy 

Afe + - 2 = Tr (Ljlk - Id ) - ^ , 

Id 

(Afc - ir + (l/Afe - If - Tr {LlL.r - 2Tr (L^L,) - , ^ 



256 



Combing with (|7.1ip we find 



^-1 UG N-1 ^6 

TJiM = npNAH] - E # ' and a% = Var[/3jvAi/] ^ ^ ^ + ^ . 

fe=i fc=i 

Moreover, the Lindenberg or simply Lyapunov condition in the general central limit theorem (see |18j ) is 
verified since we work with a sum of random variables, thus concluding the first part of the proof. 
Recalling 



we compute the convergent Riemann sums for p = 6 and p = 12. For the case ^ we have 

lim ^ y /i? = 



N^OO N " DP 

k=l 



If = we obtain 



lim -±- YhP^ hm — y 2f sinP ( —n 

k=l k=l ^ 

= StP2P j sinP {^xjdx. 



Thus for p = 6 we have that the series sums to 20St^. Then the asymptotic behavior follows from the 
assumption St]^ <^ St%, and similarly St]^N^'^ <^ St%N^ in the case = 0. □ 

Remark 7.4. The two propositions proved in this section characterize the restrictions imposed by 
the stability of the resulting scheme. In Proposition [721 stability in the large system size limit, N +oo, 
is equivalent to the inequality (|7.10p . In this case the restriction of the time-step size is imposed by the 
numerical integrator. On the other hand the stability for the scheme which uses a Metropolis corrector is 
linked to the acceptance rate of the Metropolis step. In Proposition 17.21 stability in the large system size 
limit is equivalent to the non-vanishing Metropolis acceptance rate, which is equivalent to bounded from 
above average energy variation m^r and bounded variance ctjv of the energy variation. In either case, the 
IMMP method {P > 0) induces a relative increase of order N for the boundary of numerical stability as 
compared to the exact dynamics 9 — integrated with the Verlet scheme. 

Appendix A. Surface measures. Let M"^ be endowed with the scalar product given by the positive 
definite matrix M, and consider Mz a family of sub-manifolds of co-dimension n implicitly defined by n 
independent functions Mz = {cc € | = zi, .., ^niq) = Zn} for z in a neighborhood of the origin. For 
each z in a neighborhood of the origin the conditional measure 6^i^q)^z{dq) is a measure on M.z defined in 
such a way that it satisfies the chain rule for conditional expectations with respect to the Lebesgue measure 
dq, i.e., 

dq ^ S^(^)=z{dq) dz . (A.l) 
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The surface measure ox^Mzidp) is the HausdorfF measure induced by the metric M ^ on the co-tangent 
space T*Mz = {p I V^^(g)M-V = 0}; and in the same way, cfm^ (dl) is the Hausdorff measure induced 
by the metric M on the sub- manifold Aiz- It is important to note that, although this is not explicit in the 
notation, a is defined with respect to the mass-tensor M of the mechanical system. The Liouville measure 
(JT'Mzidpdq) on the co-tangent bundle T*Mz is the volume form induced on 

T*M. = { (p, q) I ^^aq)M-'p = , aq) = 4 

by the usual symplectic form dp A dq. It can be described in terms of surface measures as follows 

CTT-Mz (dpdq) = (JT'M, (dp) CTM, (dq) . 

Finally, the co-area formula (see [T7] for a general reference) defines the relative probability density 
between S^(^q)=z{dq) and aMAd<l)- 

Proposition A.l (Co-area formula). Given the invertible Gram matrix associated with the constraints 
^(q) — z in a neighborhood of M.^ — {q \ S,{q) — 

G(g) = V^^M-iv,e, 

one has 

5^^q)=,{dq) = ^ (TAiAdq) ■ 
Vdct G{q) 

Appendix B. Langevin processes. Defining the Poisson bracket 
and the dissipation tensor 

d((7) = (Jq , 

where a is the fluctuation matrix in Definition 12.11 the Markov generator of the Langevin process in 
Definition 12.11 is 



1 



The generator £. satisfies 



C^{-,H} + ^{d,{d^,-}e-f'"}^^" 



e' 



J ipiC{ip2)e-^" dpdq^ J C*iipi)ip2e-^" dpdq, 



where 



C*={-,-H} + i{d,{d^,-}e-^"}ef^". 

The generator C* defines a Langevin process with the time-reversed Hamiltonian {—H). Reversibility of the 
process implies that the canonical measure is stationary. Furthermore, if the initial state of the system is a 
canonically distributed random variable, the probability distribution of a trajectory after the time-reversal 
is given by a Langevin process with the generator £*. When H has the form H{p, q) — ^p^AI~^p + V{q), 
reversal of impulses (p — > —p) leads to time-reversed dynamics, and a process with generator C* can be 
constructed by the following simple steps: 
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1. Reverse momenta (p -—p). 

2. Draw a random path with generator C. 

3. Reverse again momenta {p — > —p). 

When holonomic constraints, for instance, of the form 



' p^M-^Vqi =0 



e(<z) =z 

are introduced, it is useful to define the Poisson bracket on the co-tangent bundle T*M.z 

a, 6 

where T is the symplectic Gram matrix of the full constraints 

As a basic result of symplectic geometry (see [T] ) , one recovers the divergence formula with respect to the 
bracket { • , • }^ and the Liouville measure ut*m^ {dpdq) 

j {■ ^■}M,^T-MAdpdq) ^G. 

Given a constrained Langevin process in a stochastic differential equation form 

q = VpH , 

p = -VgH -jq + aW-Vg^X, 

where A are Lagrange multipliers associated with the constraints ^(g) = 0, adapted with respect to the 
noise VF, the process {pt, qt}t>o obeys hidden velocity constraints and is characterized by the stochastic 
differential equations 

q = VpH + VpS A , 

p = - VgiJ - 79 + aW - VgS A , 

where A are Lagrange multipliers associated with the full constraints q) — 0. The Markov generator 
of this process can be written in the form 

-I3h\ g/3H 



demonstrating the reversibility with respect to the constrained canonical measure e~^^ ar* M^{dp dq) . 

Appendix C. Exact sampling of fluctuation/dissipation perturbations. In this section, we 
recall how to perform exact sampling of fluctuation/dissipation perturbations. Since we only work with 
impulses, we refer to the system by using the impulse variables p only. Note that throughout the paper, we 
also use extended variables {p,Pz), however, the presentation that follows covers general cases. The kinetic 
energy of the system is ^p^Mp. We impose constraints p^ M~^Vq^ = on impulses, thus p E T*Ai and 
hence the associated orthogonal projector on T*M is 

P = Id - Vg^G-^Vj^^M-i . 
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The stochastic differential equations of motion on impulses that are integrated on a time-step interval are 

[p^M-'Wg^^Q, (Cp) 
with the usual fluctuation/dissipation relation cra"^ — 2f3~^j. The Gaussian distribution of impulses 

le-^P^'^-'PaT^,Midp) (C.2) 

is invariant under the dynamics (jC.ip . 

Proposition C.l (Exact sampling of stochastic perturbation). Given the mass matrix M, suppose 
either 6t or 7 are small enough so that the condition 

|m-i<7 (C.3) 

holds in the sense of symmetric semi-definite matrices. Let U be a centered and normalized Gaussian 
vector. Consider the mid-point Euler scheme with constraints 

St , , — 

Pn+l = Pn - y 7^'^ [Pn + Pn+l) + ^ Stall ~ Vg^ A„+i 

where Xn+i is the Lagrange multiplier associated with the constraint (Cp). The Markov kernel defined by 
the transition Pn Pn+i is reversible with respect to the Gaussian distribution (jC.2p . 

Proof. After calculating the Lagrange multiplier the expression (|C.4[) can be written as 

Pn+l = Pn - jP-/P^M-\pn + p„+i) + VStPaU . 

Consider the new variable p — [3^/'^M~^^'^p, and define the symmetric matrix 

as well as K, such that KK^ = L. In terms of these new variables we obtain from (|C.4p 

Pn+l = (Id +L)-i(Id -L)pn + 2{ld +L)-^KU. (C.5) 

Moreover, the product measure ut* Mo{dpn) <JT*Ma{dpn+i) is the measure induced on the linear subspace 
of constraints by the scalar product M^^ and the Lebesgue measure dpn dpn+i. Thus in the variables 
(PmPn+i) this measure becomes, up to a constant, the measure induced by the usual Euclidean structure. 
As a consequence the log density of the random variable {Pn,Pn+i) defined by (jC.Sp with respect to this 
latter measure is equal to 

-l\Pn\^ -l{Pn+i-iU +Ly\U ^L)p„f L-\ld +i)2(p„+i-(Id +L)-i(Id -L)p„) 



-ip^+lL"l(Id +Lfpn+1 - lplL-\ld + Lf Pn 



which is indeed symmetric between p„ and Pn+i- Hence we have shown the reversibility of the induced 
Markov kernel and consequently stationarity of the canonical Gaussian distribution. □ 
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